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Abstract 

This paper continues the work contained in two previous papers, 
devoted to the study of the dynamical system generated by a delay 
differential equation that models leukemia. Here our aim is to identify 
degenerate Hopf bifurcation points. By using an approximation of the 
center manifold, we compute the first Lyapunov coefficient for Hopf 
bifurcation points. We find by direct computation, in some zones of 
the parameter space (of biological significance), points where the first 
Lyapunov coefficient equals zero. For these we compute the second 
Lyapunov coefficient, that determines the type of the degenerate Hopf 
bifurcation. 
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1 Introduction 

We consider a delay differential equation that occurs in the study of 
periodic chronic myelogenous leukemia [11], [12] 
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Here /3q, re, 5, k, r are positive parameters. Parameter k is of the form k = 
2e _7r , with 7 also positive. We chose here to take k as an independent 
parameter, instead of 7, keeping in mind the fact that, due to its definition, 
k < 2. We do not insist here on the physical significance of the model, this 
is largely presented in [TT], [T2"j . 

As usual in the study of differential delay equations, we consider the 
Banach space B = {ip : [— r, 0] h->- R, tp is continuous on [— r, 0] } , with the sup 
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norm, and for a function x : [— r, T) \— > R, T > and a < t < T, we define 
the function x 4 € £> by Xt(s) = x(t + s). 

In [6] we proved that (pQ) with the initial condition xq = (ft, (ft E B, has 
an unique, defined on [— r, oo), bounded solution. 

The equilibrium points of the problem are 

xi = 0, X2 = {^{k-l)-l) l ' n . 
d 

The second one is acceptable from the biological point of view if and only if 

y(fc-l)-l>0, (2) 

condition that implies k > 1. 

In [6] we studied in detail the stability of equilibrium solutions of ([I]). 
The results therein are briefly recalled in Section 2. 

The second equilibrium point can present Hopf bifurcation for some 
points in the parameter space [11], [12], [TJ. In [7J we considered a typical 
Hopf bifurcation point, constructed an approximation of the center mani- 
fold in that point and found the normal form of the Hopf bifurcation (by 
computing the first Lyapunov coefficient, h). From this normal form we 
obtained the type of stability of the periodic orbit emerged by Hopf bifurca- 
tion. In Section 3 we remind the ideas concerning the approximation of the 
center manifold, necessary for the computation of l\, that were developed in 
[TJ. In order to give a visual representation of the Hopf bifurcation points 
we fix two parameters (n and (3q) and in the three dimensional parameter 
space (k, d, r) we represent the surface of Hopf bifurcation points. 

Let us denote by a the vector of parameters (/3q, n, 5, k, r). If a* is 
on the Hopf bifurcation surface and h(a*) ^ 0, the Hopf bifurcation is a 
non-degenerate one, and (x%, a*) is named codimension one Hopf point |13| . 

If at a certain a**, h(a**) = 0, we have a degenerate Hopf bifurcation 
and if we fix three of the five parameters and vary the other two parameters 
in a neighborhood of the bifurcation point in the parameter space, we obtain 
a Bautin type bifurcation [9j (if the second Lyapunov coefficient, h(ct**), 
is not zero). If h(a**) = 0, hict**) ^ then (22 >«**)> is also named a 
codimension two Hopf point [13] (since Bautin bifurcation is a codimension 
two bifurcation [9]). 

In Section 4 of the present work we present the method used by us to 
explore the existence of points with l\ = for equation (pQ). Then, for a typ- 
ical such point with l\ = 0, we develop the procedure to compute the second 
Lyapunov coefficient, For this we compute a higher order approximation 
of the center manifold. This involves solving a set of differential equations 
and algebraic computation (done in Maple). 

In Section 5 we present the results obtained for our problem by using 
the methods exposed in Section 4. We found that the considered problem 
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presents points with l\ = 0, and we identify such points (to a certain approx- 
imation) by the interval bisection technique applied with respect to one of 
the parameters. We found that all Hopf codimension two points previously 
determined have I2 < 0. 

We give tables with the values of all the parameters at the points with 
l\ = that we previously determined and we plot these points on the Hopf 
bifurcation points surface. 

Section 6 presents the conclusions of our work, while Section 7 is the 
Appendix containing the differential equations for the determination of the 
approximation of the center manifold (and their solutions). 



2 Stability of the equilibrium points 

The linearized equation around one of the equilibrium points is 

z(t) = -[B 1 + 6]z(t) + kB 1 z(t-r), (3) 
where z = x — x* , x* = x\ or x* = X2, B\ = (3'(x*)x* + (3(x*), and 

a 

P{x) = — . The nonlinear part of equation (TjQ), written in the new 

variable z, will be denoted by f(z(t),z(t — r)). The characteristic equation 
corresponding to © is 

\ + 5 + B 1 = kB 1 e- Xr . (4) 

For the equilibrium point x\ = 0, we have B\ = /3q. In [B] we pointed out 
that for xi, Re(A) < for all eigenvalues A, iff the condition fy(k — l) — l < 
is satisfied. It follows that x\ is locally asymptotically stable in this zone of 
the parameter space. 

When @f(k — 1) — 1 = 0, A = is an eigenvalue of the linearized (around 
x\) problem. For this case, by constructing a Lyapunov function, we proved 
that x\ is stable [6]. Hence, x\ is stable iff it is the only equilibrium point 
(see condition ([2]))- When the second equilibrium point occurs, x\ loses its 
stability. 

For the equilibrium point Xi-, 
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The study of stability performed in [6], that relies on the theoretical results 
of [2], reveals the following distinct situations for the stability of X2- 
I. B\ < that is equivalent to @r(k — 1) > ^r[, with two subcases: 
I. A. B\ < and 6 + B\ < 0. In this case ReX < for all eigenvalues A 
if and only if \5 + B±\ < \kB\\ and 

aiccos{{5 + Bi)/kBi) 1 

< r < — — , (61 

loo |<5 + -Bi| 
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where ujq is the solution in (0,?r/r) of the equation wcot(wr) = —(5 + B\). 

Remarks. 1. The condition 5+Bi < is equivalent to @r-(k—l)(n—k) > 
n that implies n > k, and becomes 

Po -(k-l)> U 



5 n — k 

2. The condition \6 + B\\ < \kB\\ is, in this case, equivalent to 

t' 



^§-{k — 1) > 1, that is the condition of existence of X2, hence only condition 



([6]) brings relevant information. 
LB. Bi < and 5 + Bi > 0. 

In this case, ReX < for all eigenvalues A if and only if 

f r , arccos ((5 + Bi)/kB\) ] . . 
5 + ^1 > IfcBxl or jtf + Bi < |JfeSi| and r < vv J u/ ^ j (7) 

where wo is defined as above. 

Hence, when these conditions are satisfied, x<i is locally asymptotically 
stable. 

Remarks. 1. The condition 5+B\ < is equivalent to ^9-{k—l){n—k) < 

n. 

2. Condition 5 + Bi > \kB\\ is equivalent to 

6 n ~ FPI 

If n < Fpx akwe inequality is satisfied since the expression in the 
left hand side is positive. 

If n > Fn"(> ^) ^ erj cases are possible: either ([8]) is satisfied and 
then ReX < /or every eigenvalue, or g-(k — 1) < — tj- and the stability 

n ~k+i 

condition is satisfied if r < arccos ({S+B^/kBx) ^ 
II. B 1 > that is equivalent to ^(fc - 1) < 

In this case, we can only have 6 + B\ > 0, and, by [2], .ReA < for all 
eigenvalues A if and only if kB\ < 8 + B\. But this inequality is equivalent 
to ^o^- 1 ) > 1 that is already satisfied, since x<i exists and is positive. It 
follows that if B\ > then x% is stable. 



3 Hopf bifurcation from the nontrivial equilibrium 

In the stability discussion of xi shortly presented above, in I., the case 
_ arccos((<5 + ffi)/(/c£i)) 

occurs on the frontier of the stability domain. 



(9) 
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Figure 1: Left- Surface of Hopf bifurcation points in the (6, k, r) space (n — 
2,fio = 1). Right- Projection of the surface on the plane (k, S) -figures done with 
Maple. 



Relations ([9]) and loq cot(uor) = —(5 + B{) imply 

u>o = V(^i) 2 - (<5 + £i) 2 (10) 

and that the pair Ai^ = ±ioJo represents a solution of (j3]). All others eigen- 
values have strictly negative real parts. 

Thus the points in the parameter space where relations ([9]) and (|10p hold 
are candidates for Hopf bifurcation. In order to have an image of the set of 
such points, we fixed n and /?o and in the space (k, 5, r) we represented the 
surface of equation 

arccos((<5 + £i)/(fcffi)) 

Both the numerator and the denominator of the ratio above show that the 
condition \S + Bi \ < \kB\ \ must be fulfilled. In Fig. 1 we present the image 
of the obtained surface for n = 2, /3q = 1 and a projection of the surface on 
the plane (k, S), indicating the domain of definition of the function in the 
RHS of (fTTj). 

Assuming that only one parameter varies (denote it by a) while the oth- 
ers are kept fixed, a point a* , for which (|lip is satisfied, is a nondegenerate 
Hopf bifurcation point if ^(a*) ^ and l\{a*) ^ 0, where l\ is the first Lya- 
punov coefficient, to be defined below. If l\{a*) = 0, the Hopf bifurcation is 
degenerate. 

In the sequel we study the existence of points a** in the parameter space, 
satisfying ([9]) and (fTU|) and ^i(a**) = 0, for the equilibrium point X2- In order 
to do this, we have to construct an approximation of the local invariant 
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center manifold and the approximation of the restriction of problem (JTJ) to 
the local center manifold. 



3.1 The center manifold, the restriction of the problem to 
the center manifold 

In [7J, by using the method of [3], the problem was formulated as an 
equation in the Banach space 

Bq = \ if) : [— r, 0] \- > M, ip is continuous on [— r, 0) A 3 lim ip(s) 6R>, 
I s— >o ) 

that is 

A(z t ) + d f(z t ), (12) 



dzt 
dt 

where f(ip) = f(tp(0), ip(-r)) 



d {s) 



0, fl£[-r,0), 

1, 8 = 0, 



and, for ip e B C So, 

+ 0(0)], L( 93 ) = -(J + S 1 ) v3 (0) + A: J B 1 ^(-r). (13) 

The developments of this subsection are done for a point a*, in the 
parameter space, for which relations ([9]) and (|10p are satisfied. Hence the 
problem has two eigenvalues A12 = ±iui* and all other eigenvalues have 
strictly negative eigenvalues. 

The eigenfunctions corresponding to eigenvalues Ai^ = ±ito* are given 
by ¥1,2(3) = e ±MJ s , s G [— r, 0]. Since the eigenfunctions are complex func- 
tions, we need to use the complexificates of the spaces B, Bq that we denote 
by Be, Boc, respectively. We denote by M. the subspace of Be generated 
by <pi,2(0- 

In [7], by using the ideas in [3], we constructed a projector V : Bqc ■M., 
that induces a projector defined on Be with values also in A4. The space 
Be is decomposed as a direct sum M. © ftf, where N = (J — V)Bc- We do 
not present here the details of this construction. The projection of problem 
(PQ) on A4 leads us to the equation 

du ' — ' 

— =tw*« + *i(0)/(^i« + ¥?2tZ + v), (14) 

where 

* - 1 + + ^1 - n 

UJ [l + (J + 5!)r] 2 + a;* 2 r 2 ' 1 J 

The local center manifold is, for our problem, a C°° invariant manifold, 
tangent to the space Ad at the point z = (that is x = X2), and it is the 
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graph of a function w(a*)(-) defined on a neighborhood of zero in M and 
taking values in TV. A point on the local manifold has the form 

U(fl + Uif2 + w(a*)(uipi + U(f 2 ). 

The restriction of problem (fT]) with xq = (ft to the invariant manifold is 

du — 

— = uj*iu + *i(0)/(<^iit + (p 2 u + w(a*)(u(pi + U(f2)), (16) 

with the initial condition u(0) = no, where V((ft) = UQtp\+UQip2- The real and 
the imaginary parts of this complex equation, represent the two-dimensional 
restricted to the center manifold problem. We can study this problem with 
the tools of planar dynamical systems theory (see, e.g. [9]). 



We set w(a*)(u,u) := w(a*)(u(t)(fi +u(t)(f2), and write 

w(a*)(u,u) = ^ -fn w jk(a*)u J u k , (17) 



jW 

j+k>2 J 



where Wjk(a*) G £>, u : [0, oo) h-> C. By using (|17p . we write 

f((fiu + ip 2 u + w(a*)(uipi +u<p 2 )) = ^2 ^uhk{ a *) uj7nk ■ ( 18 ) 



jlkl 

j+k>2 J 



Equation (|T6|) becomes 
du 



— =u*iu+ Y, -^93k{a*)u^u\ (19) 

j+k>2 J ' 

where 

g jk (a*) = *i(0)/ jfc (a*). (20) 

All the computation below are made at a = a* but, for simplicity, the 
parameter will not be written. 

3.2 First Lyapunov coefficient 

The first Lyapunov coefficient at a = a* is given by [9] 

h = ^2Re(ig 2 ogu +u*g 2 i). (21) 

In order to compute it, we have to compute the corresponding coefficients 
9jki °f (US])) an d, for this, we must find the coefficients of the series of / in 

We consider the function q(x) = (3(x)x, and its derivatives of order n > 2 
in X2, i.e. B n := q (jl \x2) = /3^ n \x2)x2 + nf3^ n ~ l \x2). Let us also consider 
the function F(x) = B 2 x 2 + B 3 x 3 + ... + B n x n + ... . 
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The nonlinear part of eq. (pQ) is 

f(z(t),z(t - r)) = -F(z(t)) + kF(z(t - r)), 
or, in terms of Zt, the function / of (|12j) is 

7(^) = -F(z t (0)) + A;F(^(-r)). 

We denote by <^ T(t) } the semigroup of operators generated by eq. 
I ) t>o 

(fT2j) . We have (by putting ip\ = ip, and by assuming that the initial function, 
4>, is on the invariant manifold) 

7(*) = -F([fmm + kFmtm-r)) = ( 22 ) 

= --B 2 [u(p(0) + 5^(0) + ^2o(0)n 2 + wn(0)ra+ ^ 02 (0)u 2 + ...] 2 - 
--£ 3 [^(0) +n^(0) + ^ 20 (0)u 2 + w u (0)uu + ^ 02 (0)u 2 + ...] 3 - ....+ 
+—kB 2 [ucp(-r) + uTp(-r) + -w 20 (-r)u 2 + u>n {-r)uu + -w m {—r)u 2 + -] 2 + 
+—^kB 3 [uip{-r)+ulp{-r)+-w w {-r)u 2 +w = 

= E ju f ^ K 

j+k>2 J 

The coefficients fjk involved in the expression of l\, (computed also in [7]), 
are obtained by equating the same order terms in the two series from (|22p 

f 20 = -B 2 (l-ke-^* ir ), 

fu = B 2 (k-l), 
f 02 = -B 2 (l-ke 2 "* ir ), 
/21 = -£3^20(0) - 2B 2 w n (0) + 2kB 2 e- i ^ r w ll {-r) + kB 2 e iuJ * r W2o(-r)- 

-B 3 (l-ke-^ r ), 

from where <7« will be easily obtained from (|20[) (with 'I'i(O) given in (|15p ). 
Since /21 depends on some values of W20 and W\\, we have to determine 
these functions first. 

The functions Wjk € Be are determined from the relation [13] . [10] . [5] 

^ E 7^^)^= E 7^T*^^Vi(^)+ (23) 

j+k>2 J ' ' j+fc>2 

+ E ^^ juk ^ s) + di E ^^( s V" fe > 

j+k>2 J ' ' j+k>2 J ' 
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that, by matching the same type terms yields differential equations for each 
w jk( s ), while the integration constants are determined from 

J t E -^ W M^U k + £ ^W«Vl(0)+ £ J_ g jk u^U k <p 2 (0) = 
j+k>2 J ' ' j+k>2 J ' ' j+fc>2 

j+k>2 J ' ' j+k>2 J ' ' j+k>2 J ' 

E. g. in order to determine W2o(s) we equate the terms containing u 2 
and obtain 

w 2q{ s ) = 2u*iw 2 o{s) +s>209 J (s) +g 2^p{s). 
We integrate between —r and and obtain 

^20(0) - w 20 (-r)e 2 ^ r = ^(1 - e"**) + ' (1 - e 3 -*"'). 
By equating the coefficients of u 2 in (|24p . the relation 

(2uj*i + Bi + 5) w 2 o(0) - kB 1 w 2 o(-r) = / 20 - 920 - 9o2 

results. From the two relations above, we find u?2o(0) and w 2 o{— r). 

The equations for all other functions Wij needed in the computation of 
l\ and I2 are all listed in the Appendix. They were obtained by symbolic 
computations in Maple. 

Remark. The functions Wjk, j + k = 2, were determined (for other 
problems) by many authors f]}/, in connection to the study of 

Hopf bifurcations. 

Once u)2o(— f), 1020 (0), u>n(— r), wu(0) computed, we determine the value 
of 521 and that of l\. 

The sign of the first Lyapunov coefficient l\ determines the type of the 
Hopf bifurcation [9]. That is, if l\ < 0, we have a supercritical Hopf bifur- 
cation, while if l\ > - a subcritical Hopf bifurcation. 



4 Hopf points of codimension two 

The Hopf points of codimension two are among the points of the surface 
of equation ([TT]). surface that, for n and (3q fixed can be represented in M 3 (as 
in Fig.l). Since they obey a new constraint, l\ = 0, they must lie on a curve 
on that surface. If we could write this constraint in a simple algebraical 
form, then we could attempt to represent the curve defined by and 
l\ = 0. Unfortunately, this is not possible, since the expression of l\, written 
in terms of the parameters of the problem is quite complex. 

In this situation, our strategy for finding Hopf points of codimension two 
(Bautin type bifurcation points) is the following. 
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• As in [7], we chose n* , /3q, k*, 8*, in a zone of parameters acceptable 
from biological point of view; 

• with the chosen values of the parameters we compute B^, (the value 
of B\ at these parameters) , next 

- if B^ > 0, then x 2 is stable for the chosen parameters and we stop; 

- if B\ < 0, we compute p* = 8* + B*, q* = k*B\- 

• - if \q*\ < \p*\, we stop; 

- if | q* | > \p*\, we determine to* and r* such that condition Q is 
fulfilled i.e. we set oj* = y (q*) 2 — (p*) 2 , r* = aiccos(p* / q*) /oj* ; 

• for the found values of parameters, the first Lyapunov coefficient, l\, 
is computed; 

• then we vary a parameter different of r, such that the absolute value 
of 1% decreases - we have chosen to vary 8; 

• the above computations are repeated for several values of 5, until we 
find two values S\, 82, of this parameter such that the values of l\ have 
opposite signs - if such values exist (obviously, at each computation, 
the value of r* changes, but the condition ReX* = is maintained); 

• then by using the bisection of the interval technique (with respect to 
5, starting from the interval [Si, 82]), we find a set of parameters such 
that l\ = (to a certain approximation). 

The obtained numerical results are presented in Section 5. 

4.1 Second Lyapunov coefficient 

In order to see if the Hopf points found by the above algorithm are 
points of codimension two, or have a higher degeneracy, we must compute 
the second Lyapunov coefficient in such a point. 

The formula for the second Lyapunov coefficient, at a Hopf codimension 
two point is [S] 



OJ" OJ* 2 



920931 - 5ll(4531 + 3522) - 7j502(540 + 5l3) - 530512 



1 

OJ* 3 



Re 



5211 ( 5ll(35l2 - 530) + 502(512 - Tj530) + 7^502503 J + 



5 1 

+5ii ( 502(7^530 + 3512) + ^502503 - 4 5ii530 



+ 3Im(g2o9ii)Img2i > + 
+-^4 { Im [5n5o2(520 - 3 52o5ii - 45n] + Im(g 2 ogii)[3Re(g 2 ogii) - 2|fif 2| 2 ]} • 



+ 
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We list below the expressions of the fjkS corresponding to the involved 
gjkS (not already given in Subsection 3.2): 

/so = -3B 2 w 20 (0) - B 3 + 3kB 2 e- i "* r w 20 (-r) + kB 3 e~ 3i "* r , 

fl2 = /2I) /03 = /30) 



/40 = -#2(3^2o(0) 2 + 4u; 3 o(0)) - 6B 3 w 20 (0) -B A + kB 2 {3w 20 (-r) 2 + 
+4e- iujr w 30 (-r)) + 6kB 3 e- 2tujr w 20 (-r) + kB 4 e~ 4im ', 

/31 = ~B 2 (3^1(0)^20(0) + 3w 2 i(0) + u>3o(0)) + B 3 {3w n {0) + 3w; 20 (0))- 
-B A + kB 2 (3w n {-r)w 20 (-r) + 3e~ iuiT w 2l {-r) + e iujr w 30 (-r))+ 
+kB 3 (3e~ 2iuir wu(-r) + 3w 20 (-r)) + kB 4 e~ 2u)ir , 

f 22 = -B 2 {2w n {0) 2 + 2w 12 {0) + w 02 (0)w 20 (0) + 2w 2 i(0)) + B 3 (w 02 {0)+ 
+4wn(0) + u» 20 (0)) - B 4 + kB 2 {2w n {-rf + 2e~ iuir w l2 {-r)+ 
+w 02 (-r)w 20 (-r) + 2e iuir w 21 (-r)) + kB 3 (w 02 (-r) e - 2iuir + 
+4w u (-r) + e 2wr w w {-r)) + kB 4 ; 

/04 = J 40 , hi = 1 3 \, and, finally, 

/32 = ~B 2 K 2 (0)w 3 o(0) + 6wn (0)^21 (0) + 3^22(0) + 3wi2(0)u;2o(0)+ 
+2w 31 (0)] - B 3 [Qw n (0) 2 + 3^(0) + 3^02(0)^20(0) + w 3 o(0)+ 
+6wn(0)u;2o(0) + 6r«2i(0)] + B 4 [3w 20 (0) + 6wii(0) + w 02 (0)} - 
-B 5 + kB 2 [Gwii(— r)w 2 i(— r) + 3wi 2 (—r)w 20 (-r) + 3e~ tujr w 22 (-r)+ 
+w 02 ( y -r)w 30 ( y -r) + 2e iur w 3 i{-r)} + kB 3 [3w 02 (-r)e- i ^w 2Q (-r) + 
+6e- i " jr w 11 (-r) 2 + 3e- 2 ^ r u;i2(-r) + e 2iu)r w 30 (-r)+ 
+6e iujr wn{-r)w 20 {-r) + 61021 (-r)] + kB 4 [3e^ r u> 20 (-r)+ 
+6e- iu)r w u (-r) + W Q 2 {-r)e-^ iuir ] + kB 5 e~ iuir . 

The values in and —r of the involved WijS are computed by solving the 
differential equations and by using the additional conditions listed in the 
Appendix. 

There is a function Wjk that requires a special treatment, i.e. w 2 i- As 
shown in [8], the two equations that determine w 2 i (0), w 2 i(—r) are depen- 
dent, and the value of 1021 (0) may be obtained by considering a perturbed 
problem (depending on a small e), computing the corresponding u; e 2i and 
by taking the limit as e — > 0. The formula obtained in [8] is 

, n s /2i(^i + ^2, p) - 2gu(p,w 20 ) - {g 20 + 2g 11 )(p,w n } -g 02 (p,w 2) 

W2l{0) = 2run -2rA + 2 ' 

(25) 

where (, ) is the bilinear form [3] , @J, defined on C([0,r],R)xC([-r, 0],R), 
by 
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<V>, if) = V(OMO) + kB 1 / r/>(C + r)^(C)dC, (26) 

= -2se wis , s G [-r,0] and p(() = -2Ce _w ^, C e [0,r]. 
After determining u>2i(0), we compute 1021 (—r) by using any of the two 
relations that connect the two values. All other necessary values of Wjk, 
j + k < 4 are determined by using the relations in the Appendix. 



5 Results 

As we asserted above, we explored the equilibrium point X2 with regard 
to the occurrence of Hopf degenerate bifurcations. As first step we search 
for points in the parameter space, where l\ = 0. 

We explored the zone of parameters (n, (3q, k) given by: 

(n, fa, k) e {1,1.5,2,3, ...,12} x {0.5,1,1.5,2,2.5} x {1.1, 1.2, 1.9} 

(we remind that k > 1 in order that X2 exist, and k < 2). 

For every (n*,jSg,fc*) fixed, we looked for a 5* and a r* such that 
i?eAi ) 2(a*) = and l\ = 0. 

We easily see that the case n = 1 is not interesting since, for n = 1, by 
— p (k-i)* > °> an ^ tn e points with B\ > are points where x 2 is 
stable and no Hopf bifurcations can occur. 

Interesting results were found for n = 1.5 and n = 2. For each 
Po £ {0.5, 1, 1.5, 2, 2.5}, and each k £ {1.1, 1.2, ...,1.9} we encountered 
the following behavior: at small values of 5 (i.e. 5 = 0.0001), we found 
B\ < 0, hence we could detect Hopf points. For these points l\ was proved 
to be positive, it decreased with increasing 5 and became negative for some 
value of 5. We then applied the bisection of interval technique as described 
in Section 4.1 and actually found Hopf points with l\ = 0. 

For each such points we computed the second Lyapunov coefficient and 
we found that all points with l\ = 0, determined before, have I2 < 0, hence 
they are codimension two Hopf points. 

For n = 1.5 the Hopf points occur for very large values of r (r > 60) 
and, since such values for r are not realistic (as asserted in [IT] . [12]). we 
present only the results obtained for n = 2. 

In Figures 2-6 we present the sets of parameters a** with l\{a**) = 
found by us, in tables and plotted on the surfaces of Hopf bifurcation points. 
These points were determined as the points where \ is approximately equal 
to 10" 10 . A better precision may be attained by continuing the bisection 
interval procedure further. 

By connecting the points with l\ = we determine approximately a 
curve of points with l\ = on each of the surfaces in Figs. 2-6. Let us 
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Figure 2: Hopf codimension two points for n = 2, /3q = 0.5. 




Figure 3: Hopf codimension two points for n = 2, /3q = 1. 



denote by (5 = <5(£;) the projection of this curve on the plane (k, 5). Our 
computations of l\ showed that on each of the surfaces in Figs. 2-6, l\ > 
for the points of the surface that are projected on the zone 5 < 5(k) of the 
plane (k, 5), and hence in this zone, the Hopf bifurcation is subcritical. For 
the points that are projected on the zone 5 > 5(k) we found l\ < 0, hence 
in this zone, the Hopf bifurcation is supercritical. 

Finally, for 3 < n < 12 the Hopf points found had only negative l\, and 
thus, we could not find any degenerate Hopf points. 



6 Conclusions 

The existence of Hopf codimension two points was investigated for the 
non-zero equilibrium point of the equation that models the periodic chronic 
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k 


6 


7' 


h 


1.1 


0.0137117887 


8.708438 


-0.0204 


1.2 


0.0271474053 


8.583841 


-0.0140 


1.3 


0.0403313662 


8.474054 


-0.0112 


1.4 


0.0532837222 


8.376752 


-0.0095 


1.5 


0.0660210946 


8.290117 


-0.0083 


1.6 


0.0785741932 


8.212697 


-0.0074 


1.7 


0.0909044966 


8.143320 


-0.0067 


1.8 


0.1030724022 


8.081025 


-0.0061 


1.9 


0.1150698062 


8.025013 


-0.0056 




1 08 



0.42 



Figure 4: Hopf codimension two points for n = 2, /3q = 1.5. 



k 


S 


r 


h 


1.1 


0.018282385 


6.531328 


-0.0203 


1.2 


0.036196540 


6.437880 


-0.014 


1.3 


0.053775154 


6.355540 


-0.0111 


1.4 


0.071044963 


6.282564 


-0.0093 


1.5 


0.088028126 


6.217588 


-0.0082 


1.6 


0.104743225 


6.159523 


-0.0073 


1.7 


0.121205995 


6.107490 


-0.0066 


1.8 


0.137429869 


6.060769 


-0.0060 


1.9 


0.153426408 


6.018759 


-0.0055 




Figure 5: Hopf codimension two points for n = 2, (3q = 2. 



myelogenous leukemia, presented in |12| . [11] . We searched such points in 
a zone of the 5-dimensional parameters space, of biological significance. We 
found that points a with l±(a) = actually exist and we present tables 
with values of the parameters where they occur. For each of these points we 
computed the corresponding second Lyapunov coefficient /2(a) and listed 
the obtained values in the above mentioned tables. We remark that for all 
the considered sets of parameters, I2 < 0. We also remark that the values of 
I2 are close to zero, and, at least for n = 2 and the chosen values of (3q the 
behavior of I2 as function of k is the same, i.e. I2 decreases when k increases. 
Since k = 2e _7r , it can not be greater than 2, and for k < 2 we found only 
negative values of I2 (we also considered in our computation the limit case 
k = 2 - that would imply 7 = 0- and for n = 2, k = 2 and for every /3q 
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k 


6 


r 


h 


1.1 


0.022852981 


5.225062 


-0.0203 


1.2 


0.045245675 


5.150304 


-0.0139 


1.3 


0.067218943 


5.084432 


-0.0110 


1.4 


0.088806203 


5.026051 


-0.0093 


1.5 


0.110035157 


4.974074 


-0.0081 


1.6 


0.130929032 


4.927618 


-0.0073 


1.7 


0.151507494 


4.885992 


-0.0066 


1.8 


0.171787337 


4.848615 


-0.0060 


1.9 


0.191783010 


4.815007 


-0.0055 




0.64 



Figure 6: Hopf codimension two points for n = 2, /3q = 2.5. 



considered we found still negative I2 in the points where l\ = 0). Hence 
among the considered points no higher order degeneracies can be found (in 
the zones with biological significance). 

The points where l\ = and !2 / are, if we vary two parameters, 
points of Bautin bifurcation (provided a certain non-degeneracy condition 
is satisfied [5]). A numerical investigation of the solutions of the equation 
for the parameters in and around the points with l\ = 0, is in course and 
will be soon submitted to publication. 



7 APPENDIX 

We list below the differential equations, their solutions, and the supple- 
mentary conditions for the functions w~k with j + k < 4, excepting W20 that 
is treated in Section 3. The right hand side of the differential equations as 
well as the supplementary conditions are obtained by symbolic computation 
in Maple, 
wn : 

dw n (s) _ s 

— -. —Sue +5 , n e j 

as 

wn(s) = wii(O) (e - 1) H (e -1), 

(Bi + S)wn(0) - kBiWn(-r) = f n - g n -g n . 

For Wo2, we have u>02 = ^20 . 
w 30 : 

dw 30 (s) 



ds 



3uiw 30 (s) + g3oe lws + .g 03 e lus + 3g 2 ow 2 o{s) + 3g 02 wn(s), 



w 30 (s) = ^3o(0) e 3 ^+-^g 3 o(e lws -e 3MS )+-^5o3(e-^ s -e 3 " ls )+3 52oe 3 " ls f w 20 (9)e-^ e d6+ 
2uj 4cj J 
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+3g 02 e 3ulls I* wu(6)e- 3 " i9 d6, 
Jo 

(3uii + B 1 + S)w 30 (0) - fcBiw 30 (-r) = / 30 - g 30 - g Q3 - 3g 20 w 20 (0) - 3g 02 wn(0). 
w 2 i : 
dw 21 (s) 



ds 



uiiw 21 {s)+g 21 e luls +g 12 e lu}s +2g 11 w 20 (s)+(g 20 +2g 11 )w 11 {s)+g 02 w 02 (s), 



w 21 ( S )=w 21 (0)e^ s + g 2lS e^ s + ^g 12 (e-^ s -e^ s ) + 2 9ll e^ s [' w^e^ d9+ 

Jo 

+ (.920 + 25n)e MS f w u (6)e-^ 9 d6 + g 02 e^ s f w 02 {6)e-^ e d6 , 



(u)i+B 1 +5)w 21 (0)-kB 1 w 2 i{-r) = f 2 i-g 2 i-g 12 -2g 11 w 20 (0)-(g 20 +2g 11 )wii(0)-g Q2 w 02 (0). 



w 03 : w 03 = w 30 . 
W12 : w 12 = W 2 i- 



w 40 : 

dw 40 (s) 
ds 



g 4 oe luJS +g 04 e tus +4:Ujiw4o(s)+4:g 3 ow 2 o(s)+4g 03 w 11 (s)+6g 20 w 3 o(s)+6g 02 w 21 (s), 



w 40 (s) = W 4o(0)e 4 ^ s +^4o(e^ s - e 4 ^ s ) + ^ 04 ( e - 4 " s - e 4 " is )+4 93oe 4MS f w w {e)e-^ e dO+ 
3 5 J 



+4g 03 e iuJls / Wll (e)e-^ s de+6g 20 e^ s / w 30 (e)e-^ w d6+6g 02 e^ ls / w 21 (0)e-^ i9 d6, 
Jo Jo Jo 

(iuji + Bx + 5)w 40 (0) - kBiw i0 (-r) = f 40 - g 40 -g 04 - 4g 30 w 2 o(0) - 4g 03 wn(0)- 

-6g 2 ow 30 (0) - 6g 02 w 21 (0). 

w 3 i : 

dW3 ^ = 2ujiw 31 (s) + g 31 e lus +g 13 e~ l " s +3g 21 w 20 (s) + {3g 12 + 530)^11 («)+ 
+3o3 w 02(s) + 3gnw 30 (s) + 3(g 20 + g n )w 21 (s) + 3g a2 w 12 (s), 
w 31 (s) = W3 i(0)e 2 ^+- ff31 ( e ^-e 2MS )+-^5i 3 (e-^ s -e 2 ^ s )+3.g 21 e 2MS f w 2 o{6)e~ 2uie d6+ 

UJ 6UJ J Q 

+ (.930 + 3g 12 )e 2MS f wxi{9)e~ iuie dB + ff 03 e 2 ^ s f Wo2 (0) e - 2 ^+ 
Jo Jo 

+3giie 2 ^ [* Wso(9)e- 2uie d6+(3g n +3g 20 )e 2 " is w 21 (0) e - 2 ^ 9 d9+3g 02 e 2 ^ s f w 12 (6)e- 2 " id d0 , 
Jo Jo Jo 

(2uji+B 1 +5)w 31 (0)-kB 1 w 31 (-r) = / 3 i-.93i-3i3~35 2 iw 2 o(0)-(g 30 +3g 12 )wii(0)- 
-5o3 w 02(0) - 3gnW7 3O (0) - 3(520 + ffn)^2i(0) - 33 02 wi 2 (0). 

w 2 2 : 

dW2 ^ = 92 2 e lu}s +g 22 e- luJS + 2g 12 w 2 o(s) + 2(g 21 +g 21 )w n (s) + 2g 12 w 02 (s) + 

16 



+9o2W 3 o(s) + (4gn + .g 20 ) w 2i(s) + (320 + 4g n )wi 2 (s) + g 02 w 3(s), 

W22 ( s ) = W22 - - 922 {e^ s - 1) + -g 22 (e~ ius - 1) + 2g 12 f w 2Q (9)d9+ 
to uj Jo 

+2(521+921)/ w n (6)d6 + 2g 12 [ w 02 (6)d6 + g 02 [ w 3Q (9)de+ 
Jo Jo Jo 

+ (5 20 + 4.gii) / w 2 i(0)d0+(g 2O + Ag u ) [ w 12 (9)d9 + g 02 f w 03 (9)d9, 
Jo Jo Jo 

(Bi+S)w 22 (0) - kBiw 22 (-r) = ,f 22 - g 22 +g 22 - 2g 12 w 20 (0) ~ 2(g 21 + g 21 )wn(0)- 

-2<? 12 u> 02 (0) -302^30(0) - (4gn +g 20 )w 21 (0) - (g 2Q + Ag n )w 12 (0) -g 2 w 03(0). 

W13 : W13 = w 3 i. 
W 4 : W 4 = w 40 - 
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